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Abstract 

This paper reports solution procedures for problems arising 
from the study of engine inlet wave propagation. The first 
problem is the study of sound waves radiated from cylindrical 
inlets. The second one is a quasi-one-dimensional problem to 
study the effect of nonlinearities and the third one is the study 
of nonlinearities in two dimensions. In all three problems Euler 
computations are done with a fourth-order explicit scheme. For 
the first problem results are shown in agreement with 
experimental data and for the second problem comparisons are made 
with an existing asymptotic theory. The third problem is part of 
an ongoing work and preliminary results are presented for this 
case . 
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I. Introduction 


A considerable amount of work has been done in the past to 
understand acoustic wave propagation problems. They occur in 
several situations such as sound radiation from engine inlets, 
exhausts and underwater acoustics. In this paper we are 
primarily interested in the problems arising from engine inlets, 
but the methods we develop here have the applicability to 
underwater acoustics problems too. These type of problems have 
two parts. The first one is the inlet wave propagation and the 
other is the radiation. Both have complex wave structures due to 
nonlinearities. One likes to calculate the sound pressure field 


both inside 

the 

inlet and in 

the 

atmosphere. This 

can be 

accomplished 

by 

solving both 

parts 

together, which 

becomes 


computationally complex, or by a coupling procedure by solving in 
the inlet and then in the atmosphere separately. A great deal of 
engineering literature exists on this subject, in particular for 
linear problems. This paper is intended to report a sequence of 
successes of Euler computations of both linear and nonlinear 
acoustic wave propagation from inlets. It turns out this is a 
natural way of doing these calculations, since the field 
equations are obtained from Euler equations. 

The difficulties in this approach are mostly attributable to 
treatment of boundary conditions. Particularly, the difficulties 
arise when one wants to prescribe boundary conditions 
numerically. It is known in linear wave propagation problems 
that it is difficult to prescribe farfield boundary conditions 
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numerically. The corresponding asymptotic conditions are 
Sommer f eld ' s radiation conditions which guarantees no reflections 
in the far field. In reference 2 the asymptotic behavior of the 
outgoing waves was used to extract a family of higher-order 
conditions. This holds in time domain as well as in frequency 
domain. In the frequency domain farfield conditions can be 
simulated by several methods. One such method is to combine 
integral equation solutions which is obtained in the farfield and 
solve acoustic equations in the near field. 7 An exact version of 
a similar procedure in the context of nonlocal boundary 
conditions is available in reference 9. There is another idea in 
which the numerical scheme follows the wave pattern and is called 
the wave envelope method. 1 In the time domain the work of 
reference 2 seems appropriate. A first-order method from the 
family of the radiation conditions of this work is implemented 
for the problem of sound radiation from cylindrical inlets which 
is discussed in Section III. This problem we considered in the 
linear context, but we solved the problem inside the inlet and in 
the atmosphere simultaneously. 

As far as nonlinear problems are concerned, we considered 
models of the inlet and solutions subject to the fact that the 
sound does not reflect at the open end of the inlet. The goal 
here is to examine if one can attenuate sound at the exit section 
subject to area variation of the inlet, mean flow Mach number and 
the sound source strength. A combination of these three 
parameters can be used to reduce noise from engines and these 
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facts were demonstrated in the last two decades experimentally. 
In this case we face a different nature of boundary conditions 
problem. We desire no reflection of acoustic waves at the exit 
section of the inlet. The field equations are nonlinear and the 
author is not aware of any nonlinear condition that will dictate 
no reflection for our model at the open end. The procedure for 
this class of problems discussed in Sections IV and V are due to 
linearization of the equation near the open of the inlet. In 
these linearized equations one can obtain incoming and outgoing 
characteristic variables and on the open end the incoming 
variable will be simply set to zero to obtain an approximate 
condition. Such a procedure becomes equivalent to specifying an 
impedance condition. We consider in Section IV a quasi-one- 
dimensional model and in Section V a strictly two-dimensional 
model . 

In all three classes of problems discussed here we used a 
fourth-order accurate scheme which is discussed in the next 
section. This scheme provided better results than other methods 
that we tried. In particular for the quasi-one-dimensional model 
(Section IV) , a spectral method with Chebyshev polynomials that 
we tried became very expensive. This was due to the fact that 
the Chebyshev points are concentrated near the end points of the 
inlet and as a result the calculations were inaccurate because 
the essential nonlinearities are at the center of the inlet. 
Multidomain technique is an alternative, but it is cumbersome in 
our situation to implement. Hence, we believe that the fourth- 
order method is superior for this class of problems. 
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The plan of the paper is as follows. In Section II we give a 
description of the fourth-order method and its usage in higher 
dimensions. Sections III, IV, and V contain the description of 
each problem under consideration. Finally, in Section VI we 
present some numerical calculations. 

We omit derivations of the field equations that we solve in 
our discussions and we refer readers to references 4 and 5. The 
derivation of field equations of the problem V will be reported 
in a forthcoming paper. 


II. Numerical Scheme 

As stated in the introduction the scheme used here is an 
extended version of MacCormack's method and is developed in 
reference 3. To provide a brief description let us consider a 
single equation of the form 

u fc + f x = h. (2.1) 

Let xj (1 < j < N) be equally spaced nodes in the computational 
domain. Define forward and backward flux difference operators by 

Pj< £ > ■ 7£ j - 8£ j.l + £ j±2 ' (2 ‘ 2) 

Then the scheme has four steps. From a time nAt to (n+ V 2 ) At 
it has a backward predictor and a forward corrector 
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(1) n 
u : = u . - 


, - 3 ap- (£ ", + Sh <"> 


u!J + ^ = 


u" ♦ u' 11 ♦ apt C a) + eh! 1 ’ 
J J D 3 


(2.3) 


In the next At/2 time-step it is changed to a forward predictor 
and a backward corrector stage as follows: 


u 


( 1 ) _ 


u n+ ly/ 2 

j 


+ apt f 
3 


n+ I/2 


+ $h 


n+ y /2 


n+1 


u. 


= v 2 


n+ V 2 , 
u j /z + 


U j 1) “ aP j f(1) 


+ ph 


( 1 ) 


(2.4) 


where a = At/6Ax, 3 = At/2 and the superscript (1) denotes 
predicted values. The other subscripts have the usual meaning of 
values evaluated at the indicated time step. From (2.2) one must 
notice that the fluxes are not defined at the end points, namely 
at j - 0,1 and at j = N+1, N+2 . There one can use suitable 
extrapolations. In all problems discussed here we used a third- 
order formula as follows: 


f j 4f j+l 6f j+2 + 4f j+3 f j+4 ' ~ 0 ' " 1) 


f j+l = 4f j - 6f j-l + 4f j- 2 - f j-3 ' U = N ' N+1 > 


(2.5) 


T 
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In more than one dimension we use operator splitting 
technique. Let us consider the following two-dimensional system: 

w. + F + G = H. (2.6) 

— t —x — y — 

If L x ( A t/2 ) and L y (At/2) denote symbolic solution operators to 
the one-dimensional equations 



It must be noted that the operators L x and Ly change their 
form for each At/2 interval according to (2.3) and (2.4). The 
above splitting is known to preserve second-order accuracy in 
time and does not change the spatial accuracy of the scheme. The 
stability criteria is chosen in such a way that it is common for 
both equations in (2.7) . 

The last aspect we consider here is the addition of 
artificial viscosity to resolve shock waves for the problems 
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considered in Section IV. The equations considered there are of 
conservation forms of the nature 


w fc + £(w) x = 0 . 


(2.9) 


The artificial viscous term added to (2.9) is of second-order and 
is of the form 



where v = 0(1) and p is the density. The difference form of 
this is 


5 U ' 

V — (|p I — I « - 

3x V p x 1 3x / Ax 


l p j+l ■ Oj I («j + l - Wj) 


* lp j - e 3-l 1 ( -j ' ij-l> 


This is a second-order formula. Thus we settle for less accuracy 
in the presence of shocks to obtain sharper shocks. 


Ill* Problem of Sound Radiation From Unflanged Cylindrical Inlets 
The problem studied in reference 4 is presented here. The 
same problem has been studied experimentally^ an d analytically 
using Wiener-Hopf techniques. ^ The goal here is to calculate 
radiated sound in the atmosphere subject to a given acoustic 
pressure field inside the inlet. Euler computations for this 


r 
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problem were done only in the absence of mean flow. With mean 
flow, which simulates flight situations, results are available in 
reference 7. Here we present the equations and approximate 
boundary conditions for sound radiation for an incident spin mode 
'm' (m > 0) of sound pressure wave at the left end of the inlet 
which is cylindrical (see Figure 1) . The field equations in this 
case are 


i P 

8t 


+ 


u + 
z 


+ 


v + imw 
r 


0 


3u 

TE 


+ 



3v 

TE 


+ 


ap _ 

Tr " 


0 


(3.1) 


3w im 
"3t + r~ 


P 


0 


where (u,v,w) are components of acoustic velocity in (z,r,0) 
coordinates and p is the acoustic pressure. For the plane wave 
case m = 0 and w = 0. 

The problem here then is to solve the system (2.1) subject 
to an incident field of the form 

P " e 1 * 1 " 2 J m( X m {r/a) ) eikt {3 * 2) 


together with hard wall conditions on the inlet (duct) wall and 

radiation condition in the atmosphere. In (3.2) X m are the 

zeroes of the Bessell function J' (z) and a is the radius of 

m 
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the duct and l are related to the wave number ka by 
m -* 


m 


= /(ka) 2 - A 2 . 


(3.3) 


By considering the reflected waves travelling in the -z direction 
at the left end of the inlet (3.1) and (3.2) give the following 
inflow conditions 


p + I” u t 

m 


= -2 ike 


il z / \ . , . 

m i « !lt 

IB m a 


(3.4) 


3v + \n 

at a _ v . 

m'- m 


' p 

On the wall of the inlet 


= 0, 


3P _ 


3r 


= 0 


(3.5) 


(3.6) 


which is the standard rigid wall condition. The conditions on 
the axis takes different forms depending on the incident spin 
mode number. They are 


V = 

0 

on 

r = 

0 

(m = 0) 

iw = 

0 

on 

r = 

0 

(m = 1) 

w = 

0 

on 

r = 

0 

(m> 2) . 


(3.7) 
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Fa r Field Boundary *2 



In the atmosphere the sound radiates without reflection. Thus 

one needs a nonref lective condition to be imposed at far 

distances. As we mentioned in the introduction several families 
of such boundary conditions for model problems are presented in 
reference 2. To apply such boundary conditions one needs a 
finite computational domain at distances larger than the region 
where one desires directivity measurements. See Figure 1 again 
for this purpose. Then at each point on these farfield 

boundaries of computational domain one can construct locally a 

circle of radius R from the origin (which is taken at the open 
end of the duct and on the center line) and apply radiation 
conditions of reference 2 at each point of the boundaries 
Y]_ # Y 2 and Y 3 » This condition can be shown for this 

situation to be 


- (u cos o + v sin a) + £ = 0, 


(3.8) 
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where a is the angle measured from the z-axis to a point on 
Y X , y 2 or y 3 . 

With the above formulation and with the scheme given in the 
last section the solution was started at a state of rest to 
achieve a steady time-harmonic state. Time increments are 
supplied from the inflow conditions. The solution for this 
problem will approach the solution of Helmholtz's equation times 
e ikt. (it immediate from (3.1) upon taking Fourier transform 
with respect to time that one obtains the Helmholtz equation.) 
The philosophy behind this procedure is known among 
mathematicians as the limiting amplitude principle. The 
advantage of doing this problems in time domain is that one does 
not encounter the problem of resolving or calculating the 
interior eigenvalues in the inlet. 


IV. Quasi-One-Dimensional Model 

In this section we describe a simplified nonlinear situation 
which has been popular in analyzing the nonlinear wave 
phenomena. We summarize the work appearing in references 5 and 
6. This model is depicted in Figure 2. The duct corresponds to 
the inlet situation and it has a constriction at the center. 
There is a steady flow from right to left and the sound 
propagates from a source upstream of this flow. The derivation 
of field equations used here are available in reference 5. They 
are 
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where q = 
u and p 
quantities 
components 
problem is 


*t + 

( u s q + q s u + uq) x = 0 

(4.1) 

u t + 

u u + + c^/_ 3 + 111 3— \ =0, 

s yZ s^q s T q s jjx 


Ap 

with A = A(x) in the area variation of the 

duct 

are 

acoustic velocity and density respectively. 

The 

with 

subscript s denote the steady mean 

flow 

and 

c s is the local sound speed in the flow. 

The 

then 

to solve (4.1) subject to 



u (0 , t) = f(t) 

(4.2) 


B(q,u) (l,t) = 0 

(4.3) 


Acoustic 

Source Termination 



i 1 i i i 

0 0.25 0.5 0.75 1.0 

Axial Position, x/L 


Figure 2 
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where in (4.2), f(t) denotes a source which varies in time. In 
(4.3) B is a linear operator which will give no reflections at 
the exit (i.e. r at x = 1) . As mentioned in the introduction (4.3) 
is derived upon linearization of (4.1). This process yields 


where 



The eigenvalues of this matrix are 


u„ + c and u - c . 
s s s s 


(4.4) 


The first one is positive and the other is negative. These signs 
give the characteristic directions of propagation. We form the 
matrix T from the eigenvectors so that T -1 A T is diagonal. 
Then the characteristic variables are 



(4.5) 


This gives 
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v 



u 




(4.6) 


where V]_ corresponds to the positive eigenvalue and v 2 to the 
negative one. At the right boundary Vi is the inflow variable 
which is set to zero to obtain the nonref lective boundary 
condition (4.3), i.e.. 


B (q,u) (l,t) = - _H = 0. 

q s s 

We remark here that for the fully two-dimensional problem 
discussed in the next section a similar procedure is used to 
obtain this type of boundary condition. 

As far as the source term is concerned it is driven with a 
harmonic input 

f(t) = A cos t 

where A is the amplitude of the source. 

Once again the numerical procedure described in Section II 
is applied to get accurate smooth solutions. For high Mach 
numbers and high source strengths the artificial viscosity which 
was discussed in that section was added to obtain solutions. 

As mentioned in the introduction the goal in this class of 
problems is to study the attenuation of sound pressure level at 
the exit section. It was experimentally observed that high 
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source strengths and Mach number of the mean flow are possible 
reasons. This is being demonstrated in this numerical solution 
and a discussion is available in Section VI. 

Using full Euler equations with different numerical scheme 
(lower order) is available in reference 14. Their treatment of 
boundary conditions are different from ours. 


V. Two-Dimensional Model 

The results obtained in the work described in the last 
section motivated us to seek two-dimensional effects of 
nonlinearities in duct wave propagation. This is part of our 
ongoing work. Complete details of the solution procedure will be 
reported elsewhere. Similar situations using a combination of 
asymptotic theory and numerical solutions has also been studied 
recently in reference 12. Our solution procedure is similar to 
the one in the last section. The physical configuration of this 
situation is depicted in Figure 3. 

We begin here with the field equations we use to simulate 
the nonlinear situation. These equations are derived from the 
Euler equations of the perturbed flow and then subtracted from a 
given mean state. Let u, v, p and p denote acoustic x 
component velocity, y component velocity density and pressure 
respectively. Let a subscript 's' on the quantities denote the 
corresponding mean flow. Then the field equations are given 
(5.1) . The pressure field is determined from isentropic relation 
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as given in (5.2) . 

It + TT ( p s u + u s p + pu ) + |y (p s v + v s p + P v ) = 0 

Tt ( p s u + u s p + pu J + jz ( 2 p s U s u + u s p + p s u2 + 2u s pu + ?) 

+ ly ( p s v s u + p s u s v + u s v s p + p s uv + u s pv " v s pu ) = 0 

(5.1) 



where c s = yp s /p s is the local sound speed in the flow. The 
system (5.1) has the form 


38 

3F 

3G 

Tt + 

Tx + 

3y 


(5.3) 


where 3. = (p 1# $ 2 , 0 3 ) and 


P 2 = p s u + u s p + pu 
P 3 = p s v + v s p + pv 


u / v and p are determined by 


u 



u s 0 


+ 0i 


1 


and 



P 


0 


1 * 



Figure 3 
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The area variation (see Figure 3) is included as follows. Let 
the contour of the area variation be given by 

y = d(x) . 


We introduce new variables 


5 = x 

n = 3W ' {5 ‘ 5) 

so that n = 1 will be the surface of the inlet. This change of 
variables allows us to do the computation in a rectangle, but the 
system (5.3) takes the form 

|t < d £) + (dF) + -fn (— " n - d,) = °* (5.6) 

This system is solved together with an inflow condition and a 
nonref lective condition at exit plane. On the axis the y 
component of the velocity is set to zero and on the wall normal 
component of the velocity is zero. To derive the nonref lective 
condition at the exit plane we considered the variation of (5.6) 
only in the x direction, i.e., 

ft < d l> + -f? < d I> - 


(5.7) 
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This equation is then linearized to obtain characteristic 
variables as indicated in Section IV. The nonref lect ive 
condition that comes out of this procedure is 

e 2 " (c s + u s > *1 = 0. (5.8) 

For the source plane wave incident conditions are used which 
means the pressure is prescribed. Upon linearization of (5.2) we 
have 

_ _ P(0,y,t) 
or 

B, = B±°4l±L . (5. 9) 

This completes the statement of the problem. Again the 

numerical scheme described in Section II is applied to get a 
sample solution reported in the next section. 

VI. Discussion of Results 

For the sound radiation problems discussed in Section III 
the computations were performed on a Cyber-203 machine. The 
typical grid sizes were 115 x 35. The source boundary was kept 
at eight diameters (diameter of the cylinder) away from the 
origin. The sound pressure levels were calculated on a circle at 
distance 10 diameters away from the open end. The farfield 
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boundary was chosen to enclose this circle. We compared our 
results with a result of reference 12 which was done using 
Wiener-Hopf technique. An experimental study to simulate this 
situation is available in reference 15. In this experiment a 
spinning mode synthesizer was used to produce both plane and 
spinning mode waves. A sample result comparing our results with 
both results of references 12 and 15 is presented in Figure 4. 


dB, Level 



Angle, deg 


Figure 4 


This comparison was made for a wave number ka = 3.37 and for 
the spinning mode number m = 2. In this result theoretical, ex- 
perimental and numerical results are in good agreement, except 
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near the origin for the experimental results. This is due to the 
fact that in the experiment it is difficult to completely control 
other modes and plane waves. For more comparisons we refer 
readers to reference 4. 

For the quasi-one-dimensional nonlinear model an attempt to 
compare with an experimental result (reference 8) was made. The 
comparisons were not so good. This is one of the reasons we are 
interested in the two-dimensional model. However, comparisons 
were made with the asymptotic theory of reference 10. The 
procedure discussed in Section IV was applied to a particular 
geometry called Crocco-Tsien duct. A detailed description of the 
contour of the duct is available in reference 11. This contour 
is designed in such a way that the mean flow accelerates linearly 
to Mach number one at the throat. In particular for the examples 
and comparisons given here the entry Mach number was -.50 and at 
the throat -.90 (M t ) . For this configuration the steady ohe- 
dimensional gas dynamic equations satisfied by p g and u s can be 
solved explicitly. The Euler computations are compared with the 
asymptotic theory developed in reference 10. Since the typical 
nonlinear situation arises at higher sound pressure levels and 
Mach numbers approaching one in this theory the asymptotic 
parameter was chosen on (1 - |M t |), where M t is the throat 
Mach number. A comparison is shown in Figure 5 for a throat Mach 
number -.90 and for a 140 dB source with a frequency 452 Hz. 
Solid lines are the results of our computations. The others are 
the asymptotic results at the station x/L = .75 and at 


x/L = 1.00 respectively where L is the length of the duct. 


- 22 - 


□ Asymptotic theory (ref 10), x/L = 0.75 



Shock results are also in good agreement (see reference 6) . As 
far as the sound reduction at the termination section is 
concerned a result for a higher sound pressure level (156 dB) for 
the source is presented in Figure 6. In this case acoustic 
shocks occur and cause energy loss and we can see at 5 dB sound 
pressure level drop at the termination section. 


M t = 0.90 



Figure 6 
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Around this source pressure level (156 dB) Figures 7 and 8 for 
absolute value of throat Mach numbers increasing from .70 to 90 
we have shown that the pressure level reduction is increasing and 
the wave forms starting from a smooth stage and ultimately 
becomes steeper showing the shock phenomena. Figure 7 shows 
acoustic suppression and Figure 8 shows the distortion of wave 
form. This validates the experimental suggestion that the 
increase in the Mach number attenuates sound at the termination 
or the exit section. 



Figure 7 
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Figure 8 


For the two-dimensional model we discussed in Section V we 

needed a two-dimensional flow. To simulate a situation we 

considered again the one-dimensional flow solution u„ and o 

s ►'s 

but we introduced the y component of the velocity according to 

V s (x,y) = d wr y * (6,1) 

This is valid (see reference 12) provided the variation of d(x) 
is small. For the two-dimensional case the shocks were predicted 

at low Mach numbers (reference 12) . Here we present a sample 

result that for a throat Mach number -0.85 and a 135 dB source. 
Here we observed a 10 dB sound pressure with reduction on the 
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axis at the termination section. The same pressure level 
reduction was observed in an experimental study conducted at NASA 
Langley (reference 8) . As we mentioned before the full results 
will be reported in our forthcoming paper. 


f 
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